On 
OS 



£ 



Gravitational Radiation from Coalescing Binary Neutron Stars 

Xing Zhuge, Joan M. Centrella, and Stephen L. W. McMillan 
Department of Physics and Atmospheric Science, Drexel University, Philadelphia, PA 19104 

Abstract 



q ■ We calculate the gravitational radiation produced by the merger and co- 



alescence of inspiraling binary neutron stars using 3-dimensional numerical 
simulations. The stars are modeled as polytropes and start out in the point- 



^ ■ mass limit at wide separation. The hydrodynamic integration is performed 

using smooth particle hydrodynamics (SPH) with Newtonian gravity, and the 

gravitational radiation is calculated using the quadrupole approximation. 

Q\ ■ We have run several simulations, varying both the neutron star radius 

and the equation of state. The resulting gravitational wave energy spectra 

bJQ. dE/df are rich in information about the hydrodynamics of merger and coa- 

lescence. In particular, our results demonstrate that detailed information on 

X , 

both GM/Rc 2 and the equation of state can in principle be extracted from 

the spectrum. 
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I. INTRODUCTION 

Coalescing binary neutron stars are among the most promising sources of gravitational 
waves for detection by interferometers such as LIGO and VIRGO |],|2|]. Recent studies 
|| suggest that binary inspiral due to gravitational radiation reaction, and the eventual 
coalescence of the component stars, may be detectable by these instruments at a rate of 
several per year. The inspiral phase comprises the last several thousand binary orbits and 
covers the frequency range / ~ 10-1000Hz, where the broad-band interferometers are most 
sensitive. During this stage, the separation of the stars is much larger than their radii and the 
gravitational radiation can be calculated quite accurately using post-Newtonian expansions 
in the point-mass limit ||. It is expected that the inspiral waveform will reveal the masses 
and spins of the neutron stars, as well as the orbital parameters of the binary systems [@,|5]||. 

When the binary separation is comparable to the neutron star radius, hydrodynamic 
effects become dominant and coalescence takes place within a few orbits. The coalescence 
regime probably lies at or beyond the upper end of the frequency range accessible to broad- 
band detectors, but it may be observed using specially designed narrow band interferometers 
|7|] or resonant detectors ||. Such observations may yield valuable information about neutron 
star radii, and thereby the nuclear equation of state [f^PlJTOl] ■ 

Three-dimensional numerical simulations are needed to calculate the detailed hydrody- 
namical evolution of the system during coalescence. Rather than dwell on the uncertain 
details of the physics of neutron-star interiors, most studies of this problem have opted 
simply to model the neutron stars as polytropes with equation of state 

P = Kp r = Kp 1+1 ' n , (1) 

where K is a constant that measures the specific entropy of the material and n is the 
polytropic index. A choice of n — 1 (r = 2) mimics a fairly stiff nuclear equation of 
state. Shibata, Nakamura, & Oohara |TT| , |i"2"|| have studied the behavior of binaries with 
both synchronously rotating and non-rotating stars, using an Eulerian code with gravita- 
tional radiation reaction included. Rasio & Shapiro [0,Q have simulated the coalescence 



of synchronously rotating neutron-star binaries using the Lagrangian smooth particle hydro- 
dynamics (SPH) method with purely Newtonian gravity. Recently, Davies et al. jn| have 
carried out SPH simulations of the inspiral and coalescence of nonsynchronously rotating 
neutron stars, focusing on the thermodynamics and nuclear physics of the coalescence, with 
particular application to gamma ray bursts. All of these studies use the quadrupole formula 
to calculate the gravitational radiation emitted. 

Stars in a synchronous binary rotate in the same sense as their orbital motion, with spin 
angular velocity equal to the orbital angular velocity, as seen from a non-rotating frame. 
In most close binary systems (for example, those with normal main-sequence components) 
viscosity acts to spin up initially non-rotating stars, causing them to come into a state of 



synchronous rotation in a relatively short period of time [16]. However, realistic neutron star 



viscosities are expected to be quite small, and recent work suggests |17j that the timescale 
for synchronization of neutron star binaries is generally much longer than the timescale 
for orbital decay and inspiral due to the emission of gravitational waves. Thus neutron star 
binaries are generally not expected to become synchronous as they evolve toward coalescence. 
As a complement to full 3-D hydrodynamical simulations, Lai, Rasio, & Shapiro [TE| . |T9 



have used quasi-equilibrium methods to focus on the last 10 or so orbits before the surfaces 
of the neutron stars come into contact. During this time, as tidal effects grow, the neutron 
stars are modeled as triaxial ellipsoids inspiraling on a sequence of quasi-static circular 
orbits. Using an approximate energy variational method these authors have modeled both 
synchronous and non-synchronous binaries. They find that, for sufficiently incompressible 
polytropes (n < 1.2), the system undergoes a dynamical instability which can significantly 
accelerate the secular orbital inspiral driven by radiation reaction. (This instability is driven 
by Newtonian hydrodynamics; see |2(| for the case of an unstable plunge driven by strong 



spacetime curvature.) They have calculated the evolution of binaries as they approach the 
stability limit and the orbital decay changes from secular to dynamical in character, and 
have investigated the resulting gravitational wave emission |T%|] . Their results provide an 



important component in understanding the behavior of the full 3-D hydrodynamical models. 



We have carried out 3-D simulations of the coalescence of non-rotating neutron stars 
using SPH, with particular application to the resulting gravitational wave energy spectrum 
dE/df . Our initial conditions consist of identical spherical polytropes of mass M and radius 
R on circular orbits with separations sufficiently large that tidal effects are negligible. The 
stars thus start out effectively in the point-mass regime. The gravitational field is purely 
Newtonian, with gravitational radiation calculated using the quadrupole approximation. 
To cause the stars to spiral in, we mimic the effects of gravitational radiation reaction by 
introducing a frictional term into the equations of motion to remove orbital energy at the 
rate given by the equivalent point-mass inspiral. As the neutron stars get closer together 
the tidal distortions grow and eventually dominate, and coalescence quickly follows. The 
resulting gravitational waveforms match smoothly onto the point-mass inspiral waveforms, 
facilitating analysis in the frequency domain. We focus on examining the effects of changing 
R and the polytropic index n on the gravitational wave energy spectrum dE/df. 

This paper is organized as follows. In Sec. [II] we present a brief description of the 
numerical techniques we used to do the simulations. Sec. [TTJ discusses the calculation of the 
gravitational wave quantities, including the spectrum dE/df. The use of a frictional term in 
the equations of motion to model the inspiral by gravitational radiation reaction is discussed 
in Sec. fVj and the initial conditions are given in Sec. |V|. The results of binary inspiral and 
coalescence for a standard run with M = 1.4M , R = 10km, and polytropic index n — 1 



(T = 2) are given in Sec. |V|, with the frequency analysis and the spectrum dE/df presented 
in Sec. [VI l] . Sec. |VLII] presents the results of varying the neutron star radius R and the 
polytropic index n. Finally, Sec. |IX| contains a summary of our results. 

II. NUMERICAL TECHNIQUES 



Lagrangian methods such as SPH ]2T] are especially attractive for modeling neutron-star 
coalescence since the computational resources can be concentrated where the mass is located 
instead of being spread over a grid that is mostly empty. We have used the implementation of 



SPH by Hernquist & Katz p2| known as TREESPH. In this code, the fluid is discretized into 
particles of finite extent described by a smoothing kernel. The use of variable smoothing 
lengths and individual particle timesteps makes the program adaptive in both space and 
time. 



Gravitational forces in TREESPH are calculated using a hierarchical tree method KB 



optimized for vector computers. In this method, the particles are first organized into a 
nested hierarchy of cells, and the mass multipole moments of each cell up to a fixed order, 
usually quadrupole, are calculated. To compute the gravitational acceleration, each particle 
interacts with different levels of the hierarchy in different ways. The force due to neighboring 
particles is computed by directly summing the two-body interactions. The influence of more 
distant particles is taken into account by including the multipole expansions of the cells 
which satisfy the accuracy criterion at the location of each particle. In general, the number 
of terms in the multipole expansions is small compared to the number of particles in the 
corresponding cells. This leads to a significant gain in efficiency and allows the use of larger 
numbers of particles than would be possible with methods that simply sum over all possible 
pairs of particles. 

TREESPH uses artificial viscosity to handle the shocks that develop when stars collide 
and coalesce. The code contains three choices for the artificial viscosity; we have chosen 
to use the version modified by the curl of the velocity field. This prescription reduces the 
amount of artificial viscosity used in the presence of curl, and has proved to be superior to 
the other options in tests of head-on collisions of neutron stars [fj4[] and global rotational 



instability p5 ]. Since this has already been discussed in the literature, we remark only that 



the artificial viscosity consists of two terms, one that is linear in the particle velocities (with 
user-specified coefficient a) and another that is quadratic in the velocities (with coefficient 
(3), and refer the interested reader to references |22| and [Q [see especially their equations 
(3), (5), and (6)] for details. 

As has been noted above, the neutron stars are not expected to be synchronously rotat- 
ing due to their very small physical viscosity. However in computer simulations numerical 



viscosity, either present in the method itself or added explicitly as artificial viscosity, can 
have a similar effect and cause the stars to spin up. We have monitored this effect in our 
simulations and have found it to be small; see Sec. [VT| below. 

III. CALCULATION OF GRAVITATIONAL RADIATION 

The gravitational radiation in our simulations is calculated in the quadrupole approxi- 
mation, which is valid for nearly Newtonian sources f26f . The gravitational waveforms are 
the transverse-traceless (TT) components of the metric perturbation hij, 



G_2_, 

c 4 r 



^ — 4 TT (2) 



where 

lij = P (xiXj - \5ijr 2 ) d 3 r (3) 

is the reduced (i.e. traceless) quadrupole moment of the source and a dot indicates a 
time derivative d/dt. Here spatial indices i,j = 1,2,3 and the distance to the source 
r = (x 2 + y 2 + z 2 ) 1 / 2 . In an orthonormal spherical coordinate system (r, 9, 0) with the center 
of mass of the source located at the origin, the TT part of hj has only four non-vanishing 
components. Expressed in terms of Cartesian components these are (cf. p7 ) 



he = (hx cos 2 4> + I yy sin 2 <p + I xy sin 20) cos 2 9 



+hz sin 9 — (I xz cos (f) + I yz sin 0) sin 29, 
I<H> = hx sin 2 4> + I yy cos 2 4> - I xy sin 20, (4) 

h<f> — i<i>e 

= ^(I yy - I xx ) cos 9 sin 2(p 

+I xy cos 9 cos 20 + (I xz sin — I yz cos 0) sin 9. 

The wave amplitudes for the two polarizations are then given by 



G 1 ■■ 

h+ = -7-(lee- I<p<f>), (5) 



c 4 r 



Cx 2 • ■ , . 

h x = ~j-h<t>- (6) 



c 4 r 



For an observer located on the axis at 9 = 0, 6 = these reduce to 



G 1 ■■ 

'^+,axis ~~7~\fxx yy)i \ ) 



c 4 r 



G 2 

'^x,axis ~~7~*xy \& ) 



c* r 



The angle-averaged waveforms are defined by |2^J 

(hl) = ±-Jhldn 



(hl) = ^Jhldn, (9) 



which gives 




^(hl) - 


- ±(f&) _ f(2)^2 , ±(f(2) _ x(2)\2 , J_/x(2) _ r(2)\2 , 
15 V ra zz 1 > 15 V ± yy zz 1 ' 10 V xx yy I ' 




14(1<2)\2 , 4/r(2)\2, 4/t(2)\2 
15 V J a;y ^ ~ l^xz) ~ 15 V J/« / 5 


^(hl) - 


_ lff(2) j(2)\2 , 2/r(2)\2 , 4(r(2)\2 , 4/r(2)\2 
6 V :c:r J/J// 3 v xy 1 ' 3 V xz ) ~ 3 V yz 1 ' 



(10) 



(Note that equation ([H]) corrects some typographical errors in equation (3.12) of [J27]] .) 
The standard definition of gravitational- wave luminosity is 

where there is an implied sum on z and j, the superscript (3) indicates the third time 
derivative, and the double angle brackets indicate an average over several wave periods. 
Since such averaging is not well-defined during coalescence, we simply display the unaveraged 
quantity (G/5c 5 )J^- I\a in the plots below. The energy emitted as gravitational radiation is 

AE= j Ldt. (12) 

The angular momentum lost to gravitational radiation is 



dJi __2G // (2).(3)\\ n oN 

dt ~ 5 c 5 jk \\ jm fcm // ' *■ ' 

where ej,fe is the alternating tensor. The total angular momentum carried away by the waves 



is 



AJi = f(dJi/dt)dt. (14) 

Again, we plot the unaveraged quantity (2G/5c 5 )ejjfcf^fj^ for dJi/dt. 

The energy emitted in gravitational waves per unit frequency interval dE/df is given by 



Thorne |28[ in the form 

^ = c ^(^)f(\hM)\ 2 + \hM)\ 2 ), (is) 

where r is the distance to the source and the angle brackets denote an average over all source 
angles. We define the Fourier transform h(f) of any function h(t) by 

h(f) = / h(t)e 2 ™ ft dt (16) 

J — oo 

and 

hit) = / h(f)e- 2 * ift df. (17) 



To calculate the angle-averaged quantity (|/i+| 2 + |/ix| 2 ) we first take the Fourier transforms 
of equations fl5|) and (||), to obtain 

—rh+ = I ee -I H , (18) 

P~hx = 2ig. (19) 

The Fourier transforms /i+ and /i x have the same angular dependence, given by (^), as h + 
and /i x , respectively. The angle averaging 



\h +| 2 ) = — f\h + \ 2 dn 

1 ' ' 47rJ 



|/^x| 2 ) = ^ / |/'x| :! '/0. (20) 



gives expressions analogous to (fL0|): 



— r 2 (\h l 2 \ - ±\~f {2) - f (2) l 2 4- ±\~T {2) - f 2) \ 2 4 ±\~f (2) - ~f (2) \ 2 
Q 2 ' xI'M-l / 15 1 T xx T zz\ < 15 1 W zz I "^ 10 1 «x I yy\ 



14 I jtC 2 ) 1 2 i 4_|j( 2 )|2 i 4_|r( 2 )|2 (0-\\ 



C „2/lE |2\ _ l|j( 2 ) j( 2 )|2 1 2 1 j( 2 > |2 1 4|^( 2 )|2 , 4 1 x^ 2 ) |2 



— r Z(\~h |2\ = Ilf^ _ f^|2 , 2|x^|2 , 4|jW|2 , 4, jl 
z^y^ \ I x I / gin 2/J/ I 3 1 xy I ' 3 1 xz I ' 31 yz 



We then have 

<IM/)I 2 + l^x(/)| 2 ) = <IM/)I 2 > + (|^x(/)| 2 ), (22) 

which may be substituted into expression ([T^) for dE/df. 



We use the techniques of [Q to calculate the quadrupole moment and its derivatives. 
In particular, I and I are obtained using particle positions, velocities, and accelerations 
already present in the code to produce very smooth waveforms. This yields expressions 
similar to those of Finn and Evans |29|| . However, l\j requires the derivative of the particle 
accelerations, which is taken numerically and introduces some numerical noise into L and 
dJi/dt. This noise can be removed by smoothing; see |2J] for further discussion. We have 
applied this smoothing in producing all graphs of L and dJi/dt in this paper. 

IV. MODELING INSPIRAL BY GRAVITATIONAL RADIATION REACTION 

Widely separated binary neutron stars (that is, with separation a ^> R) spiral together 
due to the effects of energy loss by gravitational radiation reaction. Once the two stars are 
close enough for tidal distortions to be significant, these effects dominate and rapid inspiral 
and coalescence ensue. In our calculations we initially place the neutron stars on (nearly) 
circular orbits with wide enough separation that tidal distortions are negligible and the stars 
are effectively in the point-mass limit. Since the gravitational field is purely Newtonian and 
does not take radiation reaction losses into account, we must explicitly include these losses 
to cause inspiral until purely hydrodynamical effects take over. 

To accomplish this, we add a frictional term to the particle acceleration equations to 



remove orbital energy at a rate given by the point-mass inspiral expression (see |15] for a 



similar approach). The gravitational wave luminosity for point-mass inspiral on circular 
orbits is |26|J 



dE 

pm ~ HI 



32 G4 ^ 3 (23) 



5 c 5 a 5 

pm 



where \i = MiM 2 /(Mi + M 2 ) is the reduced mass, M. = Mi + M 2 is the total mass of the 

system, and the subscript "pm" refers to point-mass inspiral. For equal mass stars with 

Mi = M 2 = M and separation a = £R, we find 

dE 64^WGiWy 

dt 5 G£ s \c 2 R ' l ' 

We then assume that this energy change is due to a frictional force / that is applied at 
the center of mass of each star, so that each point in the star feels the same frictional 
deceleration. Dividing the loss equally between the two stars gives 



2 dt 



(25) 

pm 



where V is the center of mass velocity of the star. Since / acts in the direction opposite to 
V this gives an acceleration 



-. f I dE 

a 



\V\ 2 

pm I y I 



(26) 



M 2M dt 
This term is added to the acceleration of every particle, so that each particle in either star 



experiences the same frictional deceleration. The net effect is that the centers of mass of the 
stars follow trajectories that approximate point-mass inspiral. This frictional term is applied 
until tidal effects dominate, leading to more rapid inspiral and coalescence; see Sec. |VI] and 
Fig. [ll] below. (Operationally, our assignment of a particle to a "star" is based simply on 
which body it happened to belong to initially. Since the frictional term is turned off before 
coalescence occurs, the question of what to do after the stars have merged does not arise.) 
The dynamics of polytropes in purely Newtonian gravity is scale free in the sense that, 
for a given polytropic index n, the results of a calculation can be scaled for any values 
of the mass M and radius R. Inspiral by gravitational radiation reaction introduces the 

10 



dimensionless parameter GM/Rc 2 , as is explicitly evident in equations fl24|) and ( pi]) for 
our frictional model of inspiral. For neutron stars, GM/Rc 2 is determined by the nuclear 
equation of state. In the calculations below, we will vary both R (and hence GM/Rc 2 ) and 
the polytropic index n. 

V. INITIAL CONDITIONS 

Since our neutron stars start out widely separated with negligible tidal interaction, they 
are modeled initially as spherical polytropes. Because the timescale for tidal effects to 
develop is much greater than the dynamical time tp for an individual star, where 



we start with stable, "cold" polytropes produced by the method discussed in |24|]. The stars 
are then placed on a circular orbit with separation a$ = £qR in the center of mass frame of 
the system in the x — y plane. Locating the centers of mass of the individual stars at (x, y) 
positions (d=oo/2, 0) initially, the stars are then given the equivalent point-mass velocities 
for a circular orbit V y = ±(M/2a ) 1 ^ 2 . 

To ensure that the stars start out on the correct point-mass inspiral trajectories, we 
also give them an initial inward radial velocity V x as follows. For point-mass inspiral the 
separation a(i) is given by [^6 



f t \ 1 / A 
a(t) = a (l - -J , (28) 

where a$ is the separation at the initial time t = and 

T <> = 4r^jk (29) 

is the inspiral time, i.e. the time needed to reach separation a = 0. For equal mass stars, 
the initial inward velocity is thus 



V = - 

x dt 



t=o 



1 da 
2~dl 

11 



= -I*!. (30) 



Since the stars have initial separation ao = £,qR this gives 

64 c fGM\ 3 

The use of the correct initial inspiral trajectory allows us to match our gravitational 
waveforms smoothly to the equivalent point-mass waveforms. This is important when ana- 



lyzing the signals in the frequency domain, as discussed in Sec. |VI]] below. 



VI. BINARY INSPIRAL AND COALESCENCE 

We take the case M = 1.4M Q and R = 10km (so GM/Rc 2 = 0.21), with polytropic 
index n — 1 and initial separation ao = 4i? as our standard model, which we refer to as 
Run 1. The parameters of this model and the other two models introduced in Sec. |VTi l 
below are presented in Table |. The results of the simulations are summarized in Table ||. 
Time is measured in units of the dynamical time to given in equation (|27|) ; for Run 1, 
t D = 7.3 x 10" 5 s. 

The evolution of this system for the case of N = 4096 particles per star is shown in Fig. |1|. 
Each frame shows the projection of all particles onto the x — y plane. As the stars spiral 
together their tidal bulges grow. By t ~ lOOtr, the stars have come into contact, after which 
they rapidly merge and coalesce into a rotating bar-like structure. Note that the merger is a 
fairly gentle process and, in contrast to head-on collisions, does not generate strong shocks 
|T4] , |24]j3~l| . Spiral arms form as mass is shed from the ends of the bar. Gravitational torques 
cause angular momentum to be transported outward and lost to the spiral arms. The arms 
expand and eventually form a disk around the central object. By t = 200to the system is 
roughly axisymmetric. 

Contour plots at t — 200t£> reveal more details of the system. In Fig. 0(a), which shows 
a cut along the x — y equatorial plane, we see that the core is essentially axisymmetric 
out to cylindrical radius w ~ 2R, where w = (x 2 + y 2 ) 1 ! 2 . As the spiral arms wind up, 
expand, and merge, the disk grows increasingly axisymmetric. In the process the arms 



12 



expand supersonically, producing shock heating that causes the disk to puff up. This can be 
seen in Fig. 0(b), which shows density contours (two per decade in density) for a cut along 
the meridional x — z plane. 

The angular velocity Q(w) of the particles is shown as a function of cylindrical radius 
w in Fig. |3|. For our choice of parameters, £1 = 1 (t^ 1 ) corresponds to a spin period 
T B pin = 0.46ms. At t = 150£d the object is in the final stage of its gravitational wave "ring 
down" (cf. Fig. || below). Fig. |||a shows that the central core w ^ 2R is still differentially 
rotating at this time (with £1 ~ m~ 0A ). The disk w > 2R is also differentially rotating, 
with £1 ~ za~ 2 . By t = 200t j o the central object has less differential rotation and is more 
nearly rigidly rotating, with £1 ~ 0.65^^, giving a spin period T spin ~ 0.71ms. The disk is 
differentially rotating with £1 ~ w^ 1 ' 7 . (Recall that Keplerian motion has £1 ~ n7~ L5 .) 

The mass m contained within cylindrical radius w is plotted in Fig. |], showing that 
~ 6% of the mass has been shed to the disk w ^ 2R. Between the times t = 150to and 
t = 200t£>, some of the matter in the disk is redistributed out to larger radii. The specific 
spin angular momentum j(-ctj) within cylindrical radius w is shown in Fig. [|. About 27% 
of the angular momentum has been shed to the disk, with continued outward transport of 
angular momentum within the disk taking place between t = lbOtp and t = 200tr>. 

The gravitational waveforms rh + and rh x for an observer on the axis at 9 = and ip = 
are shown for this run in Figs. || (a) and (b), where the solid lines give the code waveforms 
and the dashed lines the point-mass results. For the first couple of orbits after the start of 
the run (T or bit = 2Tqw) the code waveforms match the point-mass predictions. As the tidal 
bulges grow, the stars spiral in faster than they would on point-mass trajectories, leading 
to an increase in the frequency and amplitude of the gravitational waveforms (cf. [|i~8|j). 
The gravitational wave amplitudes reach a maximum during the merger of the two stars 
at t ~ 105 — llOt/), then decrease as the stars coalesce and the spiral arms expand and 
form the disk. The peak waveform amplitude (c 2 /GM)rh + ~ 0.4 corresponds to a value 
h ~ 1.4 x 10~ 21 for a source at distance r = 20Mpc (the approximate distance to the Virgo 
Cluster). By t ~ 180£d the gravitational waves have shut off and the system is essentially 
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axisymmetric. 

Fig. |7| shows (a) the gravitational wave luminosity L/L Q (where L Q = c 5 /G) , (b) the en- 
ergy AE/Mc 2 emitted as gravitational radiation, (c) dJ z /dt for the gravitational radiation, 
and (d) the total angular momentum AJ Z /J (where J is the initial total angular momen- 
tum of the system) carried by the waves. In all the gravitational-wave quantities, the code 
results (solid lines) initially track the point-mass case (dashed lines) very well, departing 
significantly from the point-mass predictions somewhat before the onset of dynamical insta- 
bility. The maximum luminosity is 1.65 x 10 _4 Lo. This may be compared with the value of 
5.3 x lCT 4 Lo found by [|4[] for the case of a head-on collision with GM/Rc 2 = 0.21; for off- 
axis collisions on parabolic orbits, the maximum obtained by those authors was 1.0 x 10~ 3 L . 
The total energy radiated away after the luminosity departs from the point-mass result by 
more than 10% is 0.032Mc 2 . Again this can be compared with 0.0025Mc 2 for a head-on 
collision and a maximum of 0.016Mc 2 for off-axis collisions obtained by [p4fl . Although the 
collisions can achieve a higher gravitational wave luminosity, they radiate less energy in the 
form of gravitational waves overall because they take place on shorter timescales than the 
inspiralling binaries. 

How sensitive are these results to the resolution of the calculation? To answer this 
question we ran the same standard model with different numbers of particles per star. Fig. |8| 
shows a comparison of the waveform rh + for the cases N = 1024, 2048, and 4096 particles 
per star. It is clear that the differences in the waveforms are small. We will see in Sec. |VI] 
below that there are only slight differences in the frequency domain as well. 

In numerical simulations viscosity, whether implicit in the numerical method or added 
explicitly as artificial viscosity, can cause problems by artificially spinning up the stars 
|T5fl . To monitor this effect in our simulations we calculated the spin angular momentum 



of each star about its center of mass and compared this with the results expected for a 
synchronously rotating star (the expected result in the limit of large viscosity). In general, 
we have found that these non-physical viscous effects always remain small in our simulations. 
For example, we ran a test case consisting of initially non-spinning stars each composed of 
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N = 1024 particles on a circular orbit of constant separation a = 4R, with artificial viscosity 
coefficients a = j3 = 0.3. After 100 tr> (~ 2.8 orbits), the stars had spin angular momenta 
< 2.3% of the synchronous value. We conclude that numerical and artificial viscosities play 
negligible roles in spinning up the stars in our simulations. 

For inspiraling stars, torquing due to the gravitationally-induced tidal bulges will cause 



a physical spin-up of the stars |1J] . This is shown for the case N = 4096 particles per star 
in Fig. |9], which plots the spin angular momentum of one star (normalized by the spin of a 
synchronously rotating star at that orbital separation) as a function of time. We see that the 
spin angular momentum of the star remains small until contact occurs at t ~ lOOtu; after 
this it increases sharply, reaching nearly 70% of the synchronous value at t = 105£/> (Each 
"star" is composed of the particles that belonged to it initially, with the orbital separation 
of the stars given by the distance between the two centers of mass.) Comparison with Fig. [I] 
confirms that this effect is due to the tidal torquing that occurs when the stars develop large 
tidal bulges, come into contact, and merge. 

Once the stars are close enough for this gravitationally-induced tidal torquing to be 

significant, Newtonian gravitational effects operating on a dynamical timescale dominate 

the secular radiation reaction effects, leading to more rapid inspiral, merger, and coalescence 

T8fl . We should turn off the gravitational friction term at some time after the Newtonian 



tidal torquing takes over and before the merger occurs, since during the merger the concept 
of equivalent point-mass trajectories is meaningless. We have experimented with turning 
off the gravitational friction term at different times and present the results for our standard 



run with N = 1024 particles per star in Fig. 10, which shows the center of mass separation 



of the two stars as a function of time. Here, the solid line shows the result of running with 
the gravitational friction term left on, and the short dashed lines show the results of turning 
this term off at (1) t = 70t£>, (2) t = 85i£>, and (3) t = lOOt^. The long-dashed line shows 
the equivalent point-mass result. 

In cases (1) and (2), the stars go into nearly circular orbits (with eccentricities appropriate 
to the inspiral radial velocity at that separation) once the frictional term is turned off. 

15 



However, the trajectory in case (3) is very similar to the result when the frictional term is 
left on, indicating that the Newtonian tidal effects are dominant by this point. The center 
of mass separation of the two stars in this case is ~ 2.5R at t — lOOtp, and then rapidly 
decreases. This result is in good agreement with the prediction of a dynamical stability limit 
at a = 2.49-R by Lai, Rasio, and Shapiro |)18|| . On the basis of these tests, we have turned 
off the gravitational friction term at t = lOOto for our standard run, and all the other plots 
in this paper for this run were done with this choice. For each of the runs reported below 
with different values of the physical parameters, we have carried out such experiments to 
determine the optimal time to turn off the gravitational friction term, since the time at 
which the Newtonian tidal effects dominate differs in each case (cf. [[18|j). 

Finally, we have experimented with the values of the artificial viscosity coefficients a and 
j3. For all runs we used the values a = j3 = 0.3 during the inspiral phase. However, since 
shocks occur during the merger, coalescence, and the expansion of the spiral arms, we ran 
some tests with different amounts of artificial viscosity during these regimes. Fig. [TTJ shows 
the waveform rh + for our standard run with N = 1024 particles per star during this phase 
for three cases: solid line, a = j3 = 0.3; short-dashed line, a = 0.3, (5 = 1; and long-dashed 
line, a = 1, (3 — 2. Not surprisingly, the amplitude of the waveform is damped as a and (3 
are increased. The low viscosity case conserves energy to ~ 2% during the period 100 — 200 
tjj (after the frictional term is turned off), indicating that the evolution of the system is not 
dominated by strong shocks. Overall, the differences in energy conservation for the three 
cases are not significant. Therefore, since the low viscosity case produces the least damping 
of the waveform, we chose to use the values a = (3 = 0.3 in all of our runs. 

VII. ANALYSIS IN THE FREQUENCY DOMAIN 

Broad-band detectors such as LIGO and VIRGO should be able to measure the gravita- 
tional waveforms of inspiraling neutron star binaries in the frequency range / ~ 10 — 1000Hz. 
Comparison of these signals with waveform templates derived from post-Newtonian analysis 
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is expected to yield the neutron-star mass M |5|U. It is important to develop techniques 
to measure the neutron-star radius R since this information, coupled with M, can constrain 



the equation of state for nuclear matter |10] 



The actual merger and coalescence stages are driven primarily by hydrodynamics and are 
expected to depend on both R and the equation of state, here parametrized by the polytropic 
index n. For most neutron-star binaries, this will take place at frequencies / > 1000Hz. In 
this regime, shot noise limits the sensitivity of the broad-band interferometers and so these 
signals may not be detectable by them 0,0|- However, a set of specially designed narrow 
band interferometers [|7| or resonant detectors || may be able to provide information about 
this high frequency region ||. 

The merger and coalescence of the neutron stars takes place within several orbits fol- 
lowing initial contact, after which the gravitational radiation shuts off fairly rapidly as the 
system settles into a roughly axisymmetric final state ||. This rapid shutoff of gravita- 
tional waves is expected to produce a sharp cutoff in the spectrum dE/df. Since the fre- 
quency of the radiation calculated in the point-mass approximation at separation a scales as 
~ a -3 / 2 ~ R~ 3 / 2 , a set of narrow-band detectors that can locate the cutoff frequency where 
the energy spectrum dE/df drops sharply may in principle determine the neutron-star radius 

#1111. 

We have calculated the spectrum dE/df for our simulations using equation (|l|). For 
point-mass inspiral, dE/df ~ / -1 / 3 p8| , where the decrease in energy with frequency reflects 



the fact that the binary spends fewer cycles near a given frequency / as it spirals in. To see 
any cutoff frequency in our data, we need a reasonably long region of point-mass inspiral in 
the frequency domain. Although our runs do start out in the point-mass regime, the binaries 
undergo dynamical instability and rapid merger within just a few orbits. To compensate 
for this we match our waveforms h + and h x onto point-mass waveforms extending back to 
much larger separations and hence lower frequencies. 

The energy spectrum dE/df for Run 1 with N = 4096 particles per star is shown in 



Fig. |T|. The solid line shows the spectrum for the extended waveform including point-mass 
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inspiral, and the short-dashed line shows the spectrum of the simulation data only. Note that 
the two curves match closely. The separation at which the data were matched corresponds 
to frequency ~ 770Hz, which is well within the inspiral regime dE/df ~ / _1//3 - Fig. |I2| 
shows that the match is smooth, and does not affect the merger and coalescence region of 
the spectrum. We have also examined the effect of using different numbers of particles on 



dE/df The result is shown in Fig. [13], which plots the spectra for Run 1 with N = 1024 
and 4096 particles per star. The use of a smaller number of particles makes only a slight 
difference to dE/df. 

Examination of Figs [12] and [13| reveals several interesting features. Starting in the point- 
mass regime, as / increases, dE/df first drops below the point-mass inspiral value, reaching 
a local minimum at / ~ 1500Hz. We identify this initial dip with the onset of dynamical 
instability. For the parameters of Run 1, Lai, Rasio, and Shapiro [|1|] found that dynamical 
instability occurs at separation a = 2A9R; for point-mass inspiral, the frequency at this 
separation is fd yn = 1566Hz. The instability causes the spectrum dE/df to drop below the 
point-mass value, since the stars fall together faster than they would had they remained on 
strictly point-mass trajectories. For reference, Fig. |T2] also shows the frequency / con tact = 
(l/2ir)(GM/ R 3 ) 1 / 2 ~ 2200Hz, which is twice the orbital frequency (in the point-mass limit) 
at separation 2R. 

This initial fall-off in the spectrum is rather slight. At higher frequencies, dE/df increases 
above the point-mass result, reaching a fairly broad maximum at / pea k ~ 2500Hz, roughly 
the frequency of the waves in Fig. |] near t ~ 125£d (the approximate time at which the 
gravitational waves shut off). To further demonstrate that this feature is associated with 
the late-time behavior of the merged system, we have calculated the spectrum dE/df for 
the cases in which the waveforms rh + and rh x (including the point-mass inspiral) were 



truncated at t — 120tn an d t = 150£d. The results are shown in Fig. [TJ], where the solid 
line shows the spectrum for the complete waveforms and the dashed lines show the spectra 
for the truncated ones. We see that this peak forms between t = 120to and t = 150£d, 
and therefore associate it with the transient, rotating bar-like structure formed immediately 



following coalescence; cf. Fig. [1} The angular speed of this structure is ~ O.GSi^ 1 (see 
Fig. ^| ), which corresponds to gravitational radiation with frequencies near ~ 2800Hz. The 
conclusion that / pea k is associated with a bar is strengthened by Run 3 below, in which the 
bar survives for a much longer time and the peak is correspondingly stronger. 

Beyond /peak, the spectrum drops sharply, eventually rising again to a secondary maxi- 



mum at / S ec ~ 3200Hz, too high to be associated with the bar. Fig. |T^ shows that this peak 
also appears between t = 120£d and t = 150£d. We attribute this broad high-frequency fea- 
ture to transient oscillations induced in the coalescing stars during the merger process — the 
result of low-order p-modes with frequencies somewhat higher than the Kepler frequency in 
the merging object (see, e.g., f33f). 

The three frequencies /d yn , /peak and / sec serve as a useful means of characterizing our 
runs. They are indicated on Fig. |12| and are presented in more detail in Table [TlJ below. 

VIII. THE EFFECTS OF CHANGING THE NEUTRON STAR RADIUS AND 

EQUATION OF STATE 

The energy spectrum dE/df shows rich structure in the frequency range / ~ 1000 — 
3000Hz in which the merger and coalescence of the neutron stars take place. To understand 
how observations of dE/df might provide information on the neutron star radius and equa- 
tion of state, we must investigate the effects of changing R and the polytropic index n. In 
this section we present the results of two runs which begin to explore this parameter space. 
We will continue this study in future papers. 

Run 2 is the same as Run 1 except that the initial neutron star radius is R = 15km. 
Taking M = 1.4M , this gives GM/Rc 2 = 0.14; see Table |. This model was run with 
iV = 1024 particles per star. The gross features of the evolution of this model are similar 
to those found in Run 1. The stars first come into contact at t ~ 250£d. By the end of the 
run, the core w ^ 2R is essentially axisymmetric and has 92% of the mass and 65% of the 
angular momentum. The disk extends out to ~ 10R. The gravitational waveforms rh + and 
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rh x are shown in Figs. [15] (a) and (b) for an observer on the axis at 9 = 0, <fi = 0. Fig. |16 
shows (a) the gravitational wave luminosity L/L and (b) the energy AE/Mc 2 emitted as 
gravitational radiation. As in Fig. [7], the time-dependence of the angular momentum carried 
away by the waves is quite similar to that of the energy, and is not presented here. In these 
figures, the solid lines give the code waveforms and the dashed lines the point-mass results. 
Some interesting properties of this model are summarized in Table [H|. 

The energy spectrum dE/df for Run 2 is shown in Fig. [17|. Again, we matched the code 
data to point-mass inspiral waveforms for analysis in the frequency domain. For Run 2 the 
match occurs at frequency ~ 420Hz. Given the parameters of this run, dynamical instability 



is expected to occur at separation a = 2A9R [0; the point-mass inspiral frequency at this 
separation is /d yn = 852Hz. Fig [TFJ shows that, as in Run 1, the spectrum drops below the 



point-mass inspiral result near / dyn . The spectrum does not then rise above the point-mass 
result at / pea k ~ 1500Hz as in Run 1; however, it does drop sharply just beyond /peak, rising 
again to a secondary peak at / sec ~ 1750Hz. See Table [TIT]. 

We estimate the frequency of the waves in Fig. |T5| at the time when the gravitational 
radiation shuts off (around t ~ 270££>) to be ~ 1300Hz — that is, close to f pea k- The orbital 
angular velocity near the end of the run is Vt ~ O.Gt^ 1 or 4.5 x 10 3 rad/s; any residual non- 
axisymmetric material rotating at this speed would yield gravitational waves at frequency 
~ 1500Hz. Again, we interpret the secondary peak as the result of high-frequency oscillations 
in the merging system. The absence of a strong peak at / pea k and the weaker maximum at 
/ sec is the result of weaker tidal forces at the point of dynamical instability, leading to a less 
pronounced and shorter-lived bar. 

Since the frequency of the gravitational radiation for point-mass inspiral is ~ a -3 / 2 ~ 
i? -3 / 2 , we expect the features in the spectrum for Run 2 to occur at lower frequencies than 
in Run 1, roughly in the ratio / Run i// Ruil 2 ~ (^Runi/^iw)" 372 ~ 1-8. Our numerical 
simulations do indeed show this behavior. For example, the ratio of the frequencies at which 
the first peak occurs is ~ 2500Hz/ 1500Hz ~ 1.7. The ratio of the frequencies at which the 
secondary peak occurs is ~ 3200Hz/1750Hz ~ 1.8. 
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Run 3 is the same as Run 1 except that we use polytropic index n = 0.5 (T = 3). This 
model was run with N = 1024 particles per star, with initial separation a = A.5R. The 
stars first make contact at t ~ lQ7to- By the end of the run, the core w ^ 2R has 93% of the 
mass and 67% of the angular momentum; the disk extends out to ~ 50R. The gravitational 
waveforms rh + and rh x are shown in Fig. [T8| for an observer on the axis at 9 = 0, = 0. 
Fig. |II| shows (a) the gravitational wave luminosity L and (b) the energy AE emitted as 
gravitational radiation. Again, the solid lines give the code waveforms and the dashed lines 
the point-mass results. Table |T| summarizes some features of this run. 

However, unlike the previous cases, the core of the merged object is slightly non- 



axisymmetric, as shown in Fig. 20. The effect of this rotating, bar-like core can be seen 



in the gravitational waveforms rh + and rh x in Fig. 18. At late times the angular velocity of 



the core is fi ~ O.St^, 1 , corresponding to a gravitational wave frequency / ~ 2200Hz. This 



agrees with the wave frequency calculated from Fig. 18 at t ~ 290t^>, confirming that the 



radiation at late times is due to the rotating core. Rasio & Shapiro |I3| also found that the 
coalescence of a synchronized binary with n = 0.5 resulted in a rotating bar-like core. 

The energy spectrum dE/df for Run 3 is shown in Fig. [Hj. Here, the match to point- 
mass inspiral waveforms occurs at frequency ~ 640Hz. Dynamical instability is expected 
to occur at separation a = 2.76R |18|| , which gives /dyn = 1342Hz. Again we see that the 
spectrum drops below the point-mass inspiral result near / dyn . The spectrum then rises to 
a sharp peak at f pea k ~ 2200Hz, drops sharply, then rises again to a secondary peak at 
/ sec ~ 2600Hz. In this model the gravitational radiation due to the rotating bar is in the 



frequency range of the first sharp peak. See Table |H. 

Lai, Rasio, and Shapiro [0 found that the onset of dynamical instability occurs at 
separation a = 2A9R for the parameters of Run 1 and at a = 2.76R for the parameters of Run 
3. From this we estimate that the ratio of frequencies at which the various spectral features 
occur should be /R un i//R un 3 ~ 1.2. Our simulations approximate this behavior. For example, 
the ratio of the frequencies at which the sharp drop occurs is ~ 2500Hz/2200Hz ~ 1.1. The 
ratio of the frequencies at which the secondary peak occurs is ~ 3200Hz/2600Hz ~ 1.2. 
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IX. SUMMARY AND DISCUSSION 

We have carried out SPH simulations of the merger and coalescence of identical non- 
rotating neutron stars modeled as polytropes. The stars start out in the point-mass regime 
and spiral together due to the effects of gravitational radiation reaction. Once the stars come 
into contact, they rapidly merge and coalesce. Spiral arms form as mass is shed from the 
ends of the central rotating bar-like structure. Angular momentum is transported outward 
by gravitational torques and lost to the spiral arms. The arms expand supersonically and 
merge, forming a shock-heated axisymmetric disk. The central rotating core becomes ax- 
isymmetric for n = 1, with the gravitational radiation shutting off rapidly after coalescence. 
For the stiffer equation of state n = 0.5, the rotating core is slightly non-axisymmetric and 
considerably longer-lived, and the gravitational waves decrease more slowly in amplitude. 

It is instructive to compare our results with other, related work. Davies et al. (1993) 
recently carried out SPH calculations very similar to ours with n ~ 0.71 (r = 2.4). Their 
results for non-rotating stars are similar to ours. Rasio & Shapiro |13| , ^4]| have performed 
SPH simulations of synchronously rotating systems. They found that polytropes with n = 1 
produce an axisymmetric core, and those with n = 0.5 yield a non-axisymmetric core, in 
agreement with our results. However, for their synchronously rotating models, the amplitude 
of the gravitational radiation drops off more rapidly after the merger than in our models. 



This effect was also seen by Shibata, Nakamura, and Oohara JXTJ] and may be due to the fact 
that the synchronously rotating stars are not spinning relative to one another when they 
merge, leading to less "ringing" of the resulting remnant. 

We have also calculated the energy emitted in gravitational waves per unit frequency 
interval dE/df. We find that the spectrum gradually drops below the point-mass inspiral 
value near the frequency at which the dynamical instability sets in; this causes the stars 
to spiral together faster than they would on point-mass trajectories. The spectrum then 
drops sharply near the frequency at which the waves from the main coalescence burst fall 
off. Finally, the spectrum rises again to a secondary peak at larger frequencies, the result of 
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oscillations that occur during the merger. 

The frequencies at which these features in the spectrum occur, as well as their ampli- 
tudes, depend on both the neutron star radius R and the equation of state specified by the 
polytropic index n. Our standard model, Run 1, has R = 10km and n — 1. When we change 
just the radius in Run 2 to R = 15km, the spectral features occur at frequencies that are 
lower by a factor ~ 1.8 and the "gentler" merger leads to a much lower amplitude in the 
energy spectrum, both near / pea k and in the secondary maximum. When instead we change 
just the polytropic index in Run 3 to n = 0.5, the features occur at frequencies that are a 
factor ~ 1.2 lower. The stiffer equation of state results in a longer-lived bar and a substan- 
tially stronger peak amplitude. Measurement of the three frequencies /d yn , /peak and / sec , 
along with the amplitudes of the spectrum there (relative to the point-mass result), thus 
may be used to obtain direct information about the physical state of the merging neutron 
stars. While the details of the peaks depend somewhat on the resolution of our simulations, 
the general results described here do not. 

The gravitational waveforms and the spectrum dE/df contain much information about 
the hydrodynamical merger and coalescence of binary neutron stars. Our results show 
that the characteristic frequencies depend on both the neutron star radii and the polytropic 
equation of state. We intend to expand our study to include the effects of both spin and non- 
equal masses, as well as gravitational radiation reaction. In particular, radiation reaction 
can be expected to affect the evolution of the rotating bar in Run 3, leading to changes in 
the spectrum dE/df. We will present the results of these studies in future papers. 
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TABLES 



Model 


R 


o-o 


GM/Rc 2 


to (ms) 


n 


r 


N 


Run 1 


10 km 


AR 


0.21 


0.073 


1 


2 


4096 


Run 2 


15 km 


AR 


0.14 


0.13 


1 


2 


1024 


Run 3 


10 km 


A.5R 


0.21 


0.073 


0.5 


3 


1024 



TABLE I. Parameters of the models are given. Both the initial radius R and polytropic index 
n have been varied. The neutron star mass is assumed to be M = 1.4M in all cases. Run 1 was 
also run with N = 1024 and N = 2048 particles, as discussed in the text. 



Model 


-'"core 


•^core 


(c 2 /GM)rh 


L/L 


AE/Mc 2 


Run 1 


94% 


73% 


0.4 


1.65 x 10~ 4 


0.032 


Run 2 


92% 


65% 


0.3 


2.12 x 10~ 5 


0.013 


Run 3 


93% 


67% 


0.4 


1.20 x 10~ 4 


> 0.042 



TABLE II. Results of the simulations are summarized. The core mass M core and angular 
momentum J cor e refer to material with cylindrical radius w <2i?. Peak values of the wave- 
form (c 2 /GM)rh and luminosity L/Lq are given; (c 2 /GM)rh ~ 0.4 corresponds to a value 
h ~ 1.4 x 10 -21 at distance r = 20Mpc. The quantity AE/Mc 2 is the energy lost to gravi- 
tational radiation after the stars depart significantly from the point-mass trajectory. It is still 
increasing at the end of Run 3 due to the rotating nonaxisymmetric core. 
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Model a dyn /dyn (Hz) / pcak (Hz) / sec (Hz) 

2500 3200 

1500 1750 

2200 2600 

TABLE III. Results of analyzing the simulations in the frequency domain are summarized. All 
frequencies refer to the spectrum dE/df; cf. Fig. O. Dynamical instability is expected to occur 



1 


2A9R 


1566 


2 


2A9R 


852 


3 


2.76R 


1342 



at separation ad yn ; these values are taken from [18|. The quantity /dyn is the gravitational wave 



frequency for point-mass inspiral at separation ad yn - 
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FIGURES 
FIG. 1. Particle positions are shown projected onto the x — y plane for Run 1 with N = 4096 

particles per star. Here, M = 1.4M , R = 10km, polytropic index n = 1, and initial separation 

ao = 4.R. The stars are orbiting in the counterclockwise direction. The vertical axis in each frame 

is y/R and the horizontal axis is x/R. 

FIG. 2. (a) Density contours are shown for a cut along the x — y plane for Run 1 with N = 4096 
particles per star at t = 200i/> The contour levels are 0.3, 0.1, 0.03, 0.01,. . . (the central density is 
0.6M/i? 3 ). (b) The same density contours as in part (a), but for a cut along the x — z plane. 

FIG. 3. The angular velocity Q(zu), where w = (x 2 + y 2 ) 1 ' 2 , is shown for Run 1 with N = 4096 
particles per star. S7 = 1 (i.e. t^ ) corresponds to a spin period T sp i n = 0.46ms. (a) t = 150i£>, (b) 
t = 200t D . 

FIG. 4. The mass fraction m(zu) is shown at t = 150£d and t = 200i£> for Run 1 with iV = 4096 
particles per star. 

FIG. 5. The specific angular momentum j(w), normalized to unity for the entire system, is 
shown at t = 150to and t = 200t£> for Run 1 with N = 4096 particles per star. 

FIG. 6. The gravitational waveforms rh + and rh x are shown for an observer on the axis at 
9 = 0, (/) = for Run 1 with N = 4096 particles per star. The solid lines give the code waveforms, 
and the dashed lines the point-mass results. 

FIG. 7. Various gravitational radiation quantities are shown for Run 1 with N = 4096 particles 
per star. The solid lines show the code results, and the dashed lines the point-mass values, (a) 
Gravitational wave luminosity L/Lq, where Lq = c 5 jG\ (b) Energy AE/Mc 2 emitted as gravita- 
tional radiation; (c) dJ z /dt; (d) The angular momentum AJ Z /J carried away by the waves, where 
J is the initial total angular momentum of the system. 
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FIG. 8. A comparison of the waveform rh + for Run 1 using N = 1024, 2048, and 4096 particles 
per star. The case N = 2048 particles per star was only run for 150i£>. 

FIG. 9. The ratio of J sp i n , the spin angular momentum of a star, to J S yno the value for a 
synchronously rotating star at that separation, is plotted versus time for Run 1 with iV = 4096 
particles per star. This graph shows that artificial viscosity does not significantly spin up the 
star during the inspiral. By t ~ lOOtu the stars are in contact, and the rapid spin up is due to 
gravitationally-induced tidal torquing during the merger process. 

FIG. 10. The separation between the centers of mass of the stars for Run 1 with N = 1024 
particles per star is shown. The solid line gives the result of leaving the gravitational friction term 
on. The short-dashed lines show the effect of turning off the friction term at (1) t = 70£/), (2) 
t = 85£d, and (3) t = 100£d. The long-dashed line shows the point-mass result. 

FIG. 11. The waveform rh + is shown for Run 1 with N = 1024 particles per star with three 
different values of the artifical viscosity coefficients a and /3. 

FIG. 12. The gravitational wave energy spectrum dE/df is shown for Run 1 with N = 4096 
particles per star. The solid line shows the result of matching to point-mass inspiral, and the 
short-dashed line is the result of using the simulation data only. The long-dashed line shows dE/df 
for point-mass inspiral. The frequency /dyn = 1566Hz is the orbital frequency (/gw = 2/ or b) a t 
which dynamical instability occurs pl|. In addition, / pea k ~ 2500Hz and / sec ~ 3200Hz mark the 
frequencies of the initial and secondary peaks, respectively. The frequency /contact corresponding 
to a circular point-mass orbit at separation 2R is also noted. 

FIG. 13. The gravitational wave energy spectrum dE/df is shown for Run 1 with N = 4096 
particles per star (solid line) and N = 1024 particles per star (dashed line). 
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FIG. 14. The gravitational wave energy spectrum dE/df is shown for Run 1 with N = 4096 
particles per star for the cases in which the waveforms were truncated at t = 120tn and t = 150tr>- 
The solid lines show the spectrum for the complete waveforms, while the dotted lines show the 
truncated cases. 

FIG. 15. The gravitational waveforms rh + and rh x are shown for an observer on the axis at 
9 = 0, 4> = for Run 2; compare with Fig. y. The solid lines give the code waveforms, and the 
dashed lines the point-mass results. 

FIG. 16. Various gravitational radiation quantities are shown for Run 2. The solid lines show 
the code results, and the dashed lines the point-mass values; compare with Fig. ^. (a) Gravitational 
wave luminosity L/Lq; (b) Energy AE/Mc 2 emitted as gravitational radiation. 

FIG. 17. The gravitational wave energy spectrum dE/df is shown for Run 2; compare with 
Fig. 0. 

FIG. 18. The gravitational waveforms rh+ and rh x are shown for an observer on the axis at 
9 = 0, (f> = for Run 3; compare with Figs. |6| and 15. The solid lines give the code waveforms, and 



the dashed lines the point-mass results. 

FIG. 19. Various gravitational radiation quantities are shown for Run 3; compare with Figs. 
and [l6|. The solid lines show the code results, and the dashed lines the point-mass values, (a) 
Gravitational wave luminosity L/Lq; (b) Energy AE/Mc 2 emitted as gravitational radiation. 

FIG. 20. Density contours are shown for a cut through the x — y plane of the central regions 
of (a) Run 1 (n = 1) and (b) Run 3 (n = 0.5). The same contour levels are plotted in both cases 
(as in Fig. 2); however, the central density in (b) is lower than the top contour in (a). 

FIG. 21. The gravitational wave energy spectrum dE/df is shown for Run 3; compare with 
Figs. ||and [17. 
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Compressed PostScript files of Figures 1-21 are available 
by anonymous ftp from zonker.drexel.edu 
in subdirectory papers/ns_coll_l . 
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